clc,clear;
close all;
%% 导入数据
X=readmatrix('data4.xlsx');
X=X(:,2:end); % 特征矩阵
Y=readmatrix('newdata.xlsx');
Y=Y(:,2);% 标签向量
futurenum=10;  %预测数
rate = mean(diff(X(end-20:end, :))); % 使用最近20年的数据计算平均增长率
% 初始化未来X矩阵
XFuture = zeros(futurenum, size(X, 2));
% 线性增长的特征
for i = [1, 2, 3, 5]  % 对第1、2、3、5列特征，假设是线性增长的
    Xpt = X(:, i);
    years = (1:size(X, 1))';
    p = polyfit(years, Xpt, 1); % 线性拟合
    fyears = (size(X, 1)+1:size(X, 1)+futurenum)';
    lgrowth = polyval(p, fyears);
    % 添加更大的随机波动
    XFuture(:, i) = lgrowth + randn(futurenum, 1) * 0.05;
end
% 平稳的特征
XFuture(:, 4) = repmat(X(end, 4), futurenum, 1) + randn(futurenum, 1) * 0.01;
%% 数据归一化
[x, psin] = mapminmax(X', 0, 1); % 将特征矩阵归一化到[0, 1]区间
[y,pcout]=mapminmax(Y',0,1);
%% 划分训练集与测试集
num = size(X, 1); % 总样本数
state=randperm(num);%随机打乱是预测具有泛化性
r = 0.8; % 训练集占比
trainnum = floor(num * r);
testnum = num - trainnum;

xtrain = x(:, state(1:trainnum))';
ytrain = y(state(1:trainnum))';

xtext = x(:, state(trainnum + 1:end))';
ytext = y(state(trainnum + 1:end))';

%% 构建SVR模型
% 训练SVR模型
Mdl = fitrsvm(xtrain, ytrain, 'KernelFunction', 'rbf', 'Solver', 'ISDA', 'KernelScale', 'auto');
%% 仿真测试
% 预测
re1 = predict(Mdl, xtrain);
re2 = predict(Mdl, xtext);
re1 = mapminmax('reverse',re1,pcout);
re2 = mapminmax('reverse',re2,pcout);
Ytrain=mapminmax('reverse',ytrain,pcout);
Ytext=mapminmax('reverse',ytext,pcout);
% 计算训练集和测试集的均方误差
mse_train = mean((re1 - Ytrain).^2);
mse_test = mean((re2 - Ytext).^2);
%% 性能评价
% 计算训练集和测试集的R^2分数
r2_train = 1 - sum((re1 - Ytrain).^2) / sum((Ytrain - mean(Ytrain)).^2);
r2_test = 1 - sum((re2 - Ytext).^2) / sum((Ytext - mean(Ytext)).^2);

%% 返回(反归一化）
re=[re1;re2];
% 计算均方根误差,R^2分数
mse_rt=sqrt(mean((re - Y(state)).^2));
r2=1 - sum((re - Y(state)).^2) / sum((Y(state) - mean(Y(state))).^2);
%% 绘图
figure;
plot(Y(state), '-*', 'LineWidth', 1, 'Color', [68, 117, 122] / 255);
hold on;
plot(re, '-o', 'LineWidth', 1, 'Color', [212, 76, 60] / 255);
legend('实际值', '预测值');
xlabel('样本');
ylabel('值');
title(['训练集预测均方误差为：', num2str(mse_rt)]);
grid;



%% 预测
newdata = XFuture; % 预测的数据
newy = newpre(newdata, psin, Mdl);
newy=mapminmax('reverse',newy,pcout);
function y = newpre(newdata, psin, Mdl)
    % 归一化
    x = mapminmax('apply', newdata', psin);
    % 预测
    y = predict(Mdl, x');
end
